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Runaway electrons, which are generated in a plasma where the induced electric field exceeds a 
certain critical value, can reach very high energies in the MeV range. For such energetic electrons, 
radiative losses will contribute significantly to the momentum space dynamics. Under certain 
conditions, due to radiative momentum losses, a non-monotonic feature - a “bump” - can form 
in the runaway electron tail, creating a potential for bump-on-tail-type instabilities to arise. Here 
we study the conditions for the existence of the bump. We derive an analytical threshold condition 
for bump appearance and give an approximate expression for the minimum energy at which the 
bump can appear. Numerical calculations are performed to support the analytical derivations. 

PACS codes: ? 


1. Introduction 


In a plasma, the drag force from Coulomb collisions acting on fast electrons decreases with 
the electron velocity. Thus, if the electric field E exceeds a threshold value E c , electrons with 
sufficient velocity will be indefinitely accelerated and are called runaway electrons. The critical 
field E c is defined as 


n e e 3 In A 
4nsbn e c 2 


( 1 . 1 ) 


where n e is the electron density, m e is the electron rest mass, c is the speed of light, e is the 
elementary charge, and In A is the Coulomb logarithm. 

Runaway electrons are generated in the presence of an induced electric field parallel to the 
magnetic field. In a tokamak, the condition E > E c can be met during the plasma start-up, during 
the flat-top phase of Ohmic plasmas if the density is sufficiently low, or in plasma disruptions. 
Especially during disruption events, a beam of runaways carrying a current of several MA and 
an energy of several MJ, may form. Such a runaway beam would pose a serious threat to the 
integrity of the first wall in reactor-size fusion devices. Any mechanism that could possibly limit 
the formation of a considerable runaway beam would be of importance. 

While the role of radiative momentum losses due to synchrotron emission has been studied 
previously in (Andersson et al. 2001), the possibility of a non-monotonic feature in the energy 
distribution of runaways - which we will henceforth refer to as a “bump” - was not considered. 
The formulation of the problem in (Andersson et al. 2001) does not ensure particle conserva¬ 
tion in the presence of radiation reaction and neglects certain terms needed to describe the bump 
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(Hazeltine & Mahajan 2004; Stahl et al. 2015). The possibility of bump formation, however, 
immediately raises the question of whether the non-monotonic behavior of the distribution could 
lead to kinetic instabilities, causing a redistribution of the runaway particles, favorable for miti¬ 
gating the potential threat to the machine. A thorough investigation of conditions favoring bump 
formation is thus needed. 

In the present paper, we use analytical calculations to investigate the runaway electron dis¬ 
tribution under the combined influence of Coulomb collisions, electric field acceleration, and 
radiative momentum losses. We show the existence of a bump and derive both a threshold con¬ 
dition for the appearance of the bump and an approximate expression for its location in parallel 
momentum space. The accuracy of the analytical estimates are then tested against numerical 
simulations carried out using CODE (Landreman el al. 2014; Stahl et al. 2015). 

The paper is organized as follows. In Sec. 2, we start by describing the particle phase-space ki¬ 
netic equation. We also discuss the transformation of the kinetic equation into the guiding-center 
phase-space and give the corresponding expression in the case of a uniform plasma. Analytical 
calculation of the condition for bump-on-tail appearance, based on the guiding-center dynamics, 
are presented in Sec. 3. A comparison of the derived conditions to numerical results is presented 
in Sec. 4, before we conclude in Sec. 5. 


2. Kinetic equation 

The kinetic equation describing the dynamics of charged particles in a plasma is 

d -t + ?r- (*/,) + ^ ■ (P fa) = ClfaJbl (2.1) 

at ox op 

where C[/ fl , //,] is the collision operator for collisions between particle species a and b, z a = (x, p) 
are the phase-space coordinates, and z a = (x, p) the equations of motion. In the Fokker-Planck 
limit, the Coulomb collision operator is given by 

ClfaJbl = — ■ ( k ab [f b ]fa - D ab [f b ] • , (2.2) 

where K a blfb\ is the friction vector and O n ;,[//,] the diffusion tensor (see Appendix A for details). 
In this paper, we do not consider contributions from large-angle collisions. 

The equations of motion for a particle with charge q and mass m combine the Hamiltonian 
motion from the electric and magnetic fields E and B, and a force F that accounts for non- 
Hamiltonian dynamics 


x =v, (2.3) 

p -qE + q v X B + F. (2.4) 


Here p = ym\ is the particle momentum, and y — 1 / V1 - v 2 /c 2 = 1 + p 2 /(mc) 2 is the rela¬ 

tivistic factor. In the case considered here, the non-Hamiltonian force is the radiation reaction 
(RR) force which was first described by Lorentz (1892) in the case of a classical non-relativistic 
point charge, and was later generalized to relativistic energies by Abraham (1905) and Dirac 
(1938). As such, the Lorentz-Abraham-Dirac (LAD) force is (Pauli 1958) 


clad 


2 2 

fr 

67T£oC 3 


3y 2 y 2 

V + _( v .y)y + -|y. 


v + l21(v-« 2 


(2.5) 


The LAD-force does however contain third order time derivatives of the particle position, which 
allows for the existence of pathological solutions. For instance, the particle velocity may grow 
exponentially in the absence of external forces (E = 0, B = 0), see e.g. (Rohrlich 2007). These 
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issues have generated discussion regarding which expression to use for the RR-force. Landau & 
Lifshitz (1975) suggested a perturbative approach in which the velocity derivatives in Eq. (2.5) 
are expressed in terms of the external force only (here the Lorentz force). Ford & O’Connell 
(1993) argue that this approach is in fact the correct one. In the paper by Spohn (2000), it is 
shown that the non-physical solutions can be avoided if the LAD-force is limited on a so-called 
critical surface, and that the resulting expressions will be equivalent to those of the perturbative 
approach. We have thus chosen to adopt the perturbative approach. Furthermore, we neglect the 
electric held in the expressions for v and v in the RR-force. This is justified since the motion of the 
particle is dominated by the magnetic held in the strongly magnetized plasmas considered here. 
An excellent discussion about the RR-force can be found in a recent review paper by Di Piazza 
etal. (2012). 


2.1. Guiding-center transformation 

Because of the v x B-terrn, the particle phase-space kinetic equation in a magnetized plasma in¬ 
cludes the rapid gyromotion time-scale which is often not interesting and is expensive to resolve 
computationally. It can, however, be eliminated using guiding-center Lie-transform perturba¬ 
tion methods. The transformation of the Hamiltonian equations of motion is one of the classical 
results in modern plasma physics (see Littlejohn 1983; Cary & Brizard 2009), and the Fokker- 
Planck collision operator has been considered in (Brizard 2004; Decker et al. 2010; Hirvijoki 
et al. 2013). The hnal step necessary to formulate our problem, the transformation of the RR- 
force, was given recently in (Hirvijoki et al. 2015). 

After the transformation, the guiding-center kinetic equation for a gyro-angle averaged distri¬ 
bution function ( F a ), including Hamiltonian motion in electromagnetic helds, the RR-force, and 
Coulomb collisions in the Fokker-Planck limit, is given by 


d{F a ) | 1 d 
dt + J dZ a 



CgcFp[(T'a)], 


( 2 . 6 ) 


where Z a form the 5D guiding-center phase-space, ff is the guiding-center phase-space Jaco¬ 
bian, Z a are the Hamiltonian guiding-center equations of motion, and ('Cj’. RR y is the contribution 
from the RR-force to the guiding-center motion. Similarly to the particle phase-space operator 
in Eq. (2.2), we can write the guiding-center Fokker-Planck collision operator in phase-space 
divergence form 


Cgcw[{F a )] - j. 




a {Fq) 
gd QZP 


(2.7) 


where (K° b gc ) and (D a f b gc ) are the guiding-center Coulomb friction and diffusion coefficients. 


2.2. Equations of motion 

We solve Eq. (2.6) in a uniform plasma, using 2D guiding-center velocity space coordinates 
Z a = (p, f), where p is the absolute value of the guiding-center momentum, <; = p\\lp is the 
pitch-angle-cosine (p\\ is the guiding-center momentum parallel to the magnetic field) and the 
guiding-center Jacobian is given by ff = p 2 . In this case, the guiding-center equations of motion 
take the simple forms 


P =qE h£ 

£=$£„(!-f 2 )/p. 


( 2 . 8 ) 

(2.9) 
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where E\\ is the electric field parallel to the magnetic field. The components of the guiding-center 
RR-force in the limit corresponding to pure synchrotron emission are (Hirvijoki et al. 2015) 


/, f p \_ ypO-f) 

V gcRRf Tr 

h* )=&-& 

V gcRR/ 


Jtr 

where the radiation reaction time-scale is defined by 

67iso(mc ) 3 


3c 


q 4 B 2 


2 ry 2 Q 2 ’ 


with r = q 2 /(Ansomc 2 ) the classical electron radius and Q = qB/(ym ) the gyro-frequency. 


( 2 . 10 ) 

( 2 . 11 ) 


( 2 . 12 ) 


2.3. Collision operator 

The particle phase-space friction and diffusion coefficients, K n b[fb\ and D a b[fbi, are expressed 
in terms of the five relativistic Braams-Karney potential functions, which are weighted integrals 
of the background distribution functions //,; (see Braams & Kamey 1989) and Appendix A for 
details. If the particle species a and b coincide, the self-collisions result in a nonlinear collision 
operator. In the present study, the particle phase-space collision operator is transformed into 
the guiding-center phase-space and linearized around a Maxwellian. The integral terms of the 
linearized collision operator are neglected and only the test particle contribution is considered. 
This choice, with some further simplifications, allows analytical solution of Eq. (2.6), which will 
be discussed in Sec. 3. 

The guiding-center friction and diffusion coefficients \K“ h J and [T>°f b gc ^ that appear in 
the guiding-center Fokker-Planck operator in Eq. (2.7) are gyro-averaged projections of their 
guiding-center push-forwarded particle phase-space counterparts. For a detailed definition of the 
guiding-center friction and diffusion coefficients, we refer to (Brizard 2004; Decker et al. 2010; 
Hirvijoki et al. 2013). The general expressions are non-trivial but, in the limit of a uniform 
plasma, the test-particle operator assuming isotropic background particle distributions becomes 
diagonal with reasonably simple non-zero components 


[Kj 

S -Vijab p. 

(2.13) 

(AU 


(2.14) 

wu 

= o f) D, t b - 

p- 

(2.15) 


The coefficients vi a b, Di m i>. and where the sub-indices / and t stand for “longitudinal” and 
“transverse” with respect to the guiding-center momentum vector, are expressed in terms of the 
five Braams-Karney potentials 'i'niu) with u = p!m a and Y a b = q 2 a q 2 b In A/(47re?) according to 


„ >n„ y(ff¥i 2«9TM 

Vl,ab = 4tt- Y ab - — -2 —— , 

nib P \ ou <r on ) 


, lf 2 y 2 d^3 | 8 rdv 4 

u du uc 2 du 

15^3 l 4 5^4 


Our guiding-center Fokker-Planck operator thus becomes 


C g cFp[<^a>] - p2 Qp 


2 1 /J7\. n d ^ Fa) 

P | Vlab P {Fa) + D Uab dp 


D 


t,ab 


p 2 dt; 




(2.16) 

'VI 00 

. 

(2.17) 


■ 

(2.18) 

t 1 " 

e i\d{F a ) 

? ’ d% ’ 

(2.19) 
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where the first term with momentum derivatives is responsible for the slowing down of fast 
particles and momentum diffusion, while the second term describes scattering in pitch-angle. 


2.4. Final expression 

For the rest of the paper, to streamline notation, we shall suppress the brackets that denote the 
gyro-averaging and the sub-index from the expression for the parallel electric field. Also, we 
will drop the particle species indices as we sum over all the background species in the collision 
operator. Thus we will have v/ = Eft v i,ab, D/ — Eft Di,ab, and A = Eft A,oZ» and our kinetic 
equation in the continuity form becomes 


dF 1 d 
dt p 2 dp 



rMW 2 ) 

Tr 



F - 


p 2 D, 


dF 

dp 



(W 2 ) 




p 2 d% j 


= 0 . 


( 2 . 20 ) 


In the following, we analyze this equation in detail. We describe the formation of a bump- 
on-tail in the electron distribution function both analytically and numerically. We also study the 
threshold conditions for the bump formation and the minimum energy of the bump location. 


3. Characteristics of a bump-on-tail feature 

The RR-force in a straight magnetic field system increases with the square of the perpendicular 
momentum, ,v 2 = s 2 ( 1 - f 2 ). As a consequence, the extent of the distribution function will, 
qualitatively, be limited in ,v ± to a region where the parallel component of the total force acting 
on an electron is positive. Electrons with higher perpendicular momenta are decelerated since the 
radiation reaction force overcomes the acceleration due to the parallel electric held. Compared 
to the case without RR-force, where the distribution function is continuously expanding in s± for 
increasing values of the parallel momentum sy = A, the limited extent of the distribution in s ± 
when the RR-force is included leads to qualitatively different dynamics. 

The width of the distribution in ,v ± is approximately constant, which means that pitch angle 
scattering is increasingly more effective at higher sy in moving the electrons to the region of 
phase space where they are decelerated. A consequence of this is that a true steady state solution 
of the kinetic equation exists and the distribution function decays exponentially in the far tail, 
something which was also observed in previous works, such as (Andersson et al. 2001). Another 
new feature is the possibility of non-monotonic behavior in the tail of the steady state distribution 
function. Note that this feature cannot be correctly described if the RR-force is not implemented 
in the phase-space divergence form that conserves the phase-space density. Therefore it is over¬ 
looked by some earlier studies. To understand the properties of the bump, and its formation, we 
will start the following analysis by assuming that a bump exists in the runaway tail, and make 
assumptions regarding its properties. These assumptions will be justified a posteriori when our 
results are compared to numerical results in Sec. 4. 

Considering a possible bump-on-tail scenario, we study Eq. (2.20) in a region where the elec¬ 
trons have high velocities compared to the electron and ion thermal speeds v » \’th,e, Vth.i- Then 
the slowing-down force is dominated by electron-electron collisions and it overshadows mo¬ 
mentum diffusion. For pitch-angle-scattering, collisions with both ions and the electron bulk are 
important. 

In the limit where the bulk populations are non-relativistic Maxwellians, we have for the fric¬ 
tion coefficient at high speeds 
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and similarly for the transverse diffusion coefficient 

1 + Z e ff n e e 4 In A 1 _ 1 + Z e ff eE c m e c 
' 2 47T£q v 2 /3 


(3.2) 


where E c is the critical electric field (Eq. 1.1), Z e g is the effective ion charge, and f3 = v/c. These 
estimates coincide with the expressions in (Andersson et al. 2001). We define the normalized 
momentum s = p/{m e c ), time r = eE c t/(m e c), radiation reaction time-scale cr _1 = eE c T r /(m e c), 
and electric field E = - E/E c , and transform the kinetic equation into a dimensionless form for 
further analysis 


OF 1 d 

- + ^r — 


dr 


ds 


s 2 \E£- 


crys( 1 




d 

+ d£ 


. .. F cr£ y 1 + Z eS dF 

s y s Ls~ di f 


= 0. 


(3.3) 


As the electric field affects only the parallel acceleration we expect the system to be strongly 
biased about £ = 1. The phase-space volume element in (,v, ^-coordinates (£J ~ s 2 ), however, 
scales nonlinearly close to £ ~ 1. A better choice for further studies close to the £ = 1 region is 
to use coordinates (sy, s ± ) which have a Jacobian ~ s ± that stays constant with respect to sy. 
The new coordinates relate to ( s , £) according to 


(3.4) 

(3.5) 

and our kinetic equation expressed with (sy, s ± ) becomes 


*11 = *£ 

*l = * yjl -£ 2 - 


9F +E 9F 2F yl ( S|1 dF + S± dF \ 
dr <9i'y s s 2 \ s <9iy s ds ± ) 


y( 1 + Zeff) 

1 9 ( s 

id 2 F 

d 2 F' 

„ *11' *-l d 2 F 

7 ^ 

[ dF H 

h dF )] 

2s 

*± ds ± \ ± ds ± ) s 2 

^ll 

ds l. 

S 2 <9sy<9.Sj_ 

S 2 

\ ds w 

*|| ds ± l 


cr 

7 


(2+4ii)F + s ± (l + ii)|^ 

OS\ 


*11*1 


dF 
ds ii 


= 0. 


(3.6) 


Instead of attempting to solve Eq. (3.6), in the following we will concentrate on the dynamics at 
jj_ = 0, which will be sufficient to prove the existence of a bump and to estimate its location in 
the electron tail. 

We assume the distribution to be a smooth function of s ± , which allows us to create a power 
series expansion around s ± = 0: 


F (*II>*l) = Y 

n =0 


„2n 


(2 n)\ 


g(2n)p 

“ „2h 

, J V ^ 

g(2n+l)p 

ds ( l n) 

2j(2)! + 1)! 

(i||,0) n=0 v ’ 

ds (2n+1) 


(3.7) 


J (J||.0) 


Because the electric field is acting only in the parallel direction F is ’’even” in i.e., we can 
formally state that fT.vy, ,y ± ) = f’(,vy. -s ± ) although our phase-space does not extend to ,y ± < 0. 
Thus, all the odd s ± -derivatives at ,y ± = 0 must vanish, and we find 


F (*I1>*l) = Y 

n =0 


„2 n 


(2 n)\ 


g(2n)p 

ds ( l n) 


(i||.0) 


(3.8) 
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With the help of the expansion, we may accurately calculate the limit 


lim ——— 
s±-»o os j_ 


/ dF\ ^ 

\d 2 F 


ds]_ 


(3.9) 


Oll-O) 


and write the = 0 limit of the kinetic equation as 

(1+Zeff)^ 


dF 


dr 


(%0) 


1 


l + si 


E - 


dF 


Sin 


(3.10) 


(*[|,0) 


(1+Zeff)^ 


1 + S 2 


*11 


d 2 F 

2 

/ 

cr 

1 

ds\ 

(S||,0) 

, V 1+ 

+ ^11 

/ 


[^]( i |i,0) - 0. 


Assuming that a steady state solution exists, the possible extrema are characterized by the con¬ 
dition 


dF 


ds\\ 


= 0. 


(3.11) 


(J||,0) 


We thus find an algebraic equation that defines the locations of these extrema 


cr 

\ 

1 

(l + Z eff ) ^ 

jl + sj 

1 d 2 F 

, V 1+ 

1 s\\ 

*n 

F ds]_ 


= 0. 


(3.12) 


(«n ,ot 


3.1. Threshold condition for the appearance of the bump 

Considering a steady-state solution to Eq. (3.10), in a situation where the bump is on the verge of 
appearing, a single inflection point exists in the distribution function instead of local maxima or 
minima. In this section we derive a threshold condition describing the appearance of an inflection 
point, by requiring the first and second sy-derivatives of the distribution to vanish simultaneously. 

Before we start the analysis, we note that the steady state distribution function represented by 
Eq. (12) of (Andersson el al. 2001) is separable in ,V| and s±, and it is of the formoc exp[-W^.s 2 /2], 
where W = 2 cr/(E - 1). To find this result they neglect f/s\\ corrections compared to df /ds\\ 
terms in the kinetic equation, which is appropriate in the very far tail [sy corresponds to p\\ in the 
notation of (Andersson el al. 2001)]. Therefore, the quantity 


W 2 (sy) = - 


1 d 2 F 
F ds\ 




(3.13) 


should approach W 2 . in the ,vy —» oo limit. Thus it is useful to define /c(.s’y), so that W 2 (s\\) = 
k(s\\)W 2 , and k —> 1 as sy —> oo. Numerical calculations tell us that 0 < k < 1 for the regions of 
interest in the runaway tail, and it is often slowly varying function of ,vy. That is, the characteristic 
width of the distribution function in the s ± direction, l/W 2 , decreases with increasing ,vy, and 
slowly asymptotes to a constant value. For now, we may simply use 0 < k < 1 as a working 
hypothesis to be verified through numerical calculations later. 

We start with the Eq. (3.12) satisfied at extrema or inflection of the distribution function and 
rewrite it as 


L(s\\) = 2 (crsy + ^1 + s 2 ) - AXsyXl + s 2 ) = 0, 


(3.14) 
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where /f(sy) = W^,(l + ZeffVdVy) = E~ 1 ctk(s\\) with E = (E - 1 )/[2( 1 + Z e g)]. It is useful to form 

L'(s\\)/2 ~ cr + *!_ -K(su)s h (3.15) 

Vw 

where prime denotes a derivative with respect to sy, and we neglected a term, -K'(s\\)( 1 + sjj)/2, 
assuming k to be sufficiently slowly varying function of sy. From Eq. (3.15) we see that L'(sy) —> 
-oo as sy —» oo, while L'[ S||= o = 2cr > 0. It can also be shown that L' = 0 only has one root for 
positive values of sy, which then has to correspond to a single maximum of L. 

If the distribution has an inflection point, both L and L' should vanish there. Assuming the 
slowly varying a- to be a constant, the system of equations (3.14) and (3.15) can be solved for sy 
and K to find 


and 


K 0 = K(s w o) =(4V2 o-T 


8<r, 4- 


8 


3 + Vl + 8cr 2 


(l + 4cr 2 - Vl + 8cr 2 ) 


8 


3 + Vl + 8<x 2 


(3.16) 


%=Jl - r -r, (3-17) 

\ 3+ Vl + 8cr 2 

where the subscript 0 refers to values of quantities at the threshold of the bump appearance. By 
inspecting the expressions (3.16) and (3.17) we find that both of them increase with cr mono- 
tonically; K 0 between 2 and oo, and Syo between 0 and 1. This means that an inflection point is 
always located below ,vy = 1. Note, that we assume the inflection point to be sufficiently far from 
the bulk, and that k' is small; violation of these assumptions may move ,vyo above unity. However, 
since this problem cannot be addressed until the more complete Eq. (2.20) is solved, we assume 
that the conditions above are fulfilled in order to proceed analytically. Since sy < 1, we can make 
use of the expansion 

s 2 

^T+^=l + -l+0 ( «y). (3-18) 

By neglecting <9(.s-jj) terms (which is reasonably good even when Sy approaches unity), Eq. (3.14) 
becomes quadratic in ,vy: 

2(cr.vy + 1 + sjj/2) - AXsyXl + sj|) = 0. (3.19) 

At the inflection point syo, Eq. (3.19) must have a single root, which requires the discriminant to 
vanish. This determines the threshold value of K 

K 0 = (3 + Vl +4cr 2 ) /2, (3.20) 

which is positive. This can be substituted back into Eq. (3.19) to find 

•syo = cr/(K 0 - 1) = tr[(l + Vl + 4 cr 2 ) /2]"‘ . (3.21) 

Combining /C(^y) = E^ 1 ctk(s\\) with Eq. (3.20) to solve for a positive cr that corresponds to K > 2, 
yields cr as a function of E at the threshold for bump formation 

3 k/E + V8 +k 2 /E 2 
2 [k 2 /E 2 - l) 


O’ o = 


(3.22) 
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where k = K(sy = syo) < 1 is treated as a parameter. When cr is increased above cro, L becomes 
negative, and no bump appears. Reducing k below unity increases the threshold value of cr. Thus 
Eq. (3.22) with k = 1 represents an absolute lower threshold in cr for a monotonic behavior of 
the steady state distribution function. Furthermore, even at k = 1 the threshold is limited to the 
region E < 1, thus a bump should always appear when E > 1. 

We have considered K(s\\ — 0) > 2, in which case L(sy) = 0 can have zero (no bump), one 
(threshold), or two (bump exists) positive real roots. When ^(sy = 0) < 2 there is only a single 
positive root of L(sy). Since the distribution function cannot have positive slope in the high ,V| 
limit this root should also correspond to a bump. That, however, requires the existence of a 
minimum in the distribution function along the positive s\\ axis, which must then appear outside 
the domain of validity of Eq. (3.12). In fact, this minimum will appear close to the bulk part of 
the electron distribution, where neglected corrections to the collision operator become important. 

We can conclude that for E > 1, there should always be a bump in the steady state distribution 
function as long as there is a finite magnetic field. This, perhaps somewhat counter-intuitive, 
result needs some clarification to accommodate the well known cr = 0 limit. When no loss 
mechanisms are considered (in this case when cr = 0), the electron distribution has no steady 
state solution, and the runaway tail at = 0 should converge to a 1/sy decay. When cr is small, 
the bump location moves to high values in sy, as will be shown in the next section. The runaway 
tail always builds up starting from the bulk and, for a tiny cr, the process may take such a long 
time that the distribution never becomes non-monotonic in practice. In this scenario the steady 
state distribution and the bump have no relevance. Also, other loss mechanisms may limit the 
distribution function to momenta below syo in realistic cases. 


3.2. An estimate for the bump location in the far-tail 

In order to proceed and estimate the location of the bump and the shape of the distribution func¬ 
tion, we look for a steady-state solution in a region where the guiding-center parallel momentum 
is large. Using the expansion 


yi+Sj^iy+OCSy 1 ), 


(3.23) 


Eq. (3.10) gives 


(£- 1 +6>(V)) 


dF 


dsn 


-(1 +Zeff) 


(*11,0) 


d 2 F 


ds]_ 


~ [F\s h o) = 0 . 

(*||,0) II 


(3.24) 


We neglect the 0( Sy 1 )-term, which is valid if E -1 is not very small, and assume that the width of 
the distribution function in the ,v ± -direction, and thus W 2 = -E ffj 1 F/ds 2 ± )\( s , 0 ), stays approx¬ 
imately constant close to the bump. This essentially means that we are looking for a separable 
solution of the form F ~ hisfgis^). We obtain an ordinary differential equation 


[E - l) S\\W - 2(1 + cr)h = -(1 + Z eff )W 2 s {l h, (3.25) 

which is solved by 

h(s l{ ) ~ sJ (1+tr) /( £ -i) exp [-W 2 (l + Z dT )/(E - l)jy], (3.26) 

and the location of the bump-on-tail is given by 
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Figure 1. Typical examples of non-monotonic runaway electron distribution functions, a) The pitch angle 
average of the distribution function with (solid curve) and without (dashed) synchrotron radiation reac¬ 
tion. A Maxwellian distribution is also indicated (dash-dotted), b) Contour plot of the distribution function 
corresponding to the solid curve in (a), as a function of jy and 


If we again assume that k does not exceed unity, recalling Wy = 2 cr/(E - 1), we find a lower 
bound for the parallel momentum at the bump 


1 + cr E — 1 


1 + 0 ", 


■^11.min 


-2 E. 


(3.28) 


cr 1 + Z eff 

We see that for small values of cr, the bump would appear at high parallel momenta. By setting 
sy.min to some upper limit of physical interest, sy,, Eq. (3.28) may be used to find an estimate for 
a lower “practical limit” in cr for the appearance of the bump. Namely, if cr is smaller than 

1 


o" L 


( %l k)/( 2£) - r 


(3.29) 


for k = 1 , then a bump would only appear at some large parallel momentum above sp,, which is 
then deemed physically irrelevant. Note that if the bump is in the far tail, k can be significantly 
less than unity, as will be shown in the next section, using numerical simulations. Letting k < 1 
increases the practical limit in cr. Another implication of Eq. (3.29) is that for a normalized 
electric field higher than E — S||,l/2, the bump always appears above sy L for any value of cr. 


4. Comparison to numerical results 

The numerical results shown in this section were performed with the continuum simulation 
tool, CODE, used in its time-independent mode. CODE solves the two dimensional momentum 
space kinetic equation in a homogeneous plasma, using a linearized Fokker-Planck operator valid 
for arbitrary electron energies. For a detailed description of the tool, see (Landreman et al. 2014). 

First, we provide a typical example of a non-monotonic runaway distribution function in the 
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presence of radiation reaction. Figure la) shows the momentum dependence of the pitch angle 
averaged runaway electron distribution with ( B = 2.5 T) and without (B - 0 T) radiation reaction 
force, plotted with solid and dashed curves respectively. Technically, the pitch-angle-averaged 
distribution is the lowest mode in a Legendre polynomial expansion of F in i : . normalized so that 
F is unity at its maximum. The simulations were performed with the parameters T— 5 keV, n e — 
2 • 10 19 m \ Z e ff = 1.2, and E = 2. Note that the distribution function without radiation reaction 
represents a quasi-steady state. The lack of loss mechanisms leads to a slow but steady depletion 
of the bulk electron population, as more and more electrons run away and leave the computational 
domain. This outflow must be balanced by an artificial source of thermal (Maxwellian) electrons 
to maintain the quasi-steady state. In the presence of radiation reaction, the distribution is a 
true steady state. When the radiation reaction is included, the non-monotonic feature is present 
when the distribution is averaged over pitch angles in the present example. However, we note 
that for less pronounced bumps, the pitch-angle-averaged distribution can have a monotonic tail, 
or it may exhibit a bump at some s which are appreciably lower than those where the bump is 
observed in the full 2D-distribution. This may have an impact on the possibility of bump-on-tail 
type instabilities to arise. 

Figure lb) shows a contour plot in sy-Sj. momentum space of the distribution function cor¬ 
responding to the solid curve in Fig. la). Although this example is representative of a typical 
runaway electron distribution, the location and the height of the bump, and the width of the dis¬ 
tribution in Sj_ can vary significantly depending on the plasma parameters. The relation between 
the location sy and the local “width” [~ 1/W 2 = 1 /(Wy.K)] of the distribution given by Eq. (3.27) 
is accurate as long as the location of the bump is not close to unity, i.e. sufficiently far from 
the no-bump threshold (3.22). This justifies the approximations applied to the collision operator 
in our analysis. In particular, energy diffusion can be neglected since no sharp features of the 
distribution function in .vy are present, as seen in Figure lb). 

In order to investigate the validity of our analytical calculations, we have performed a nu¬ 
merical analysis of the appearance of the bump by scanning the parameter space with CODE. 
The electron temperature and density where held constant at the values T e = 1 keV and n e = 
5 • 10 18 m 3 , respectively, while the magnetic field, the induced electric field and the effective 
ion charge were varied over the ranges Be [1,6] T, E e [2,14] and Z e g e [1,3]. The numerical 
calculations used 950 momentum grid points, 130 Legendre modes for the decomposition in 
and a highest resolved momentum of s = 34, provided well converged solutions. 

The results of the scan are presented in Fig. 2, where circles and crosses correspond to distribu¬ 
tions with and without a bump, respectively. The color coding of the circles reflects the location 
of the bump, with 100% in the color bar corresponding to sy = 34. Simulations with a bump 
appearing above 80% (sy = 27) are excluded from the figure, since those results may be affected 
by the bump being too close to the highest resolved momentum. As expected from Eq. (3.28), 
increasing E or decreasing cr moves the bump towards larger momenta. 

A reasonably good agreement is found between the numerical calculations and the analytical 
threshold for the bump to exist. The solid curve shows this threshold, Eq. (3.22), for k = 1. 
Above cr ~ 0.5 the “no bump” solutions obey the analytical threshold and fall to the left of 
the threshold curve. There are some solutions with bump in this region as well, however this is 
not surprising, since k at the bump is allowed to be less than unity, in which case the threshold 
moves towards lower values of E. The threshold begins to fail for lower values of cr, showing 
that the approximation K' <s {cr, ,vy) used in Eq. (3. 15) breaks down. Nevertheless, the qualitative 
behavior of the threshold is still captured by Eq. (3.22). The lower right corner of the plot (high 
E and small cr) is not populated, since some simulations where the bump would have appeared at 
too high i’y were excluded. With sy ^ = 27, a k value as low as 0.3 is needed in order for the limit 
given by Eq. (3.29) (dashed line) to correspond well to the boundary of the region of excluded 
points. Thus, k can be significantly lower than unity at a bump with large momentum. 
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Figure 2. Parameter scan of CODE simulations yielding steady state solutions with (colored circles) or 
without (black crosses) a bump in the runaway tail. Good correlation is found with the analytical threshold 
condition given in Eq. (3.22) for k = 1 (solid line). The dashed line represents the “practical threshold” in 
Eq. (3.29) for sp. = 27 and k = 0.3. The color coding shows the location of the bump relative to Sy = 34. 



2E(l + cr)/a 


Figure 3. Parallel momentum of the bump. Circles denote the locations of the bumps according to numerical 
solutions, while the solid line represents a theoretical lower limit, Eq. (3.28). In the simulations, all bumps 
appear above the analytical threshold condition. 


From the parameter scan used to generate Fig. 2, the simulations exhibiting bumps were com¬ 
pared to the theoretical lower bound for the location of the bump (Eq. 3.28). As shown in Fig. 3, 
we find that, indeed, the parallel momenta at the bumps (shown with green circles) are all higher 
than the lower bound (solid line). In fact, most of the sy values are well above this limit. This 
merely confirms our finding that k at the bump is typically less than unity, especially for param¬ 
eters sufficiently exceeding the no-bump threshold. 











5. Conclusions 
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We have analyzed the runaway electron distribution function, accounting for the radiation reac¬ 
tion force, and have shown that the steady state runaway distribution can become non-monotonic. 
Furthermore, a threshold condition for the appearance of the bump, as well as a lower limit to 
its location in momentum space, were derived. While slowing down and pitch-angle scattering 
due to Coulomb collisions are taken into account in our analysis, we do not consider the effect 
of large angle collisions, and we restrict the study to a straight magnetic field geometry. Our 
analytical results show good agreement with numerical simulations obtained using the CODE 
solver. 

We find that for a normalized electric field larger than unity, E > 1, the steady state electron 
distribution always exhibits a bump, independently of the value of the cr parameter quantifying 
the strength of the radiation reaction, as long as the magnetic field is non-vanishing (cr > 0). For 
a smaller electric field, the appearance of the bump is well correlated with the cr threshold given 
by Eq. (3.22). Although above this threshold there must always be a bump in the steady state 
distribution function, it may not have a practical relevance in some cases. When cr is small and/or 
E is large, the bump would be located at very large parallel momentum, but the forefront of the 
electron distribution can require a long time to reach that far. This motivates the introduction of 
another “practical” threshold condition, Eq. (3.29). If cr is lower than this threshold, the bump will 
appear at a momentum above some specified limit, sy l, and can then be considered unimportant. 
In particular, above a normalized electric field of E = S||,l/2, this criterion is satisfied for any cr. 

Nevertheless, when the radiation reaction is strong enough and/or the parallel electric field is 
not too high, there is a possibility for a bump to form in the runaway tail. This non-monotonic 
feature presents a potential source for bump-on-tail instabilities, which can play a role in limiting 
the formation of large runaway beams. 

The authors are grateful to M. Landreman and P. Helander for fruitful discussions. IP was 
supported by the International Postdoc grant of Vetenskapsradet. 


Appendix A. The relativistic collision operator 

The relativistic particle phase-space Beliaev-Budker Collision operator in the Landau form is 
defined as 


rff n Fab 8 f i 'ut ir d f b f d f a \ 

CIA.AI = . j rfp U(»,u) -[f„ w - A ¥ ). 

where F,,/, = e 2 e 2 In A/(4jre^), u = p /m a , u' = p '/nib, and the collision kernel is given by 


(Al) 


with the coefficients 


— - w 2 1 - uu - u'u' + r(uu' + u'u) , 

L J 

(A 2) 

y — Vl + (u/c) 2 , 

(A3) 

y' =sj 1 + (m'/c) 2 , 

(A 4) 

r -yy' - u ■ u'/c 2 . 

(A 5) 

w =c Vr 2 - 1. 

(A 6) 


Braams & Karney (1989) found a corresponding differential form for the collision operator 


d 


df a 


C[fa,fb\ = -AT- • I KAMfa - D ab [f„] • ^ , 


(A 7) 
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where K a b[fb\ is the friction vector and TSablfb] the diffusion tensor that are defined with differ¬ 
ential operations on Braams-Karney potentials T„(m) according to 


Ka b [f b \ = -4tt^^K|T 1 

m y 


T 2 


= -4tt^ 

r 

The differential operators K and L are defined as 

<9T 




/ uu\ <3“T r uu\ 1 t uu\( 


<9u ’ 
<3 2 T 


dT 

du 


and the potential functions are given by the integrals 

1 

4 n , 


'Fq(u) = 




fdu'*™, 

J y'w 

f 

J y'w 


^2(U) 


1 r,, . 

= - — d u si 
8 n J 


sinh l {w/c) 


cfb( uQ 
Y 


T3(u) = -to 


^ 4 (U) = 


1 

yin 


r (? 

J du — [r sintT 1 (w/c) - (w/c))/*(u')- 


Furthermore, the potential functions satisfy differential relations 

Lo'Vo = fb, 

L\ x ¥\ =f b , 

Ltf 2 = ^i, 

T2T3 = T () , 

T 2 T 4 = t 3 , 

where the operator L* is defined 

<9 2 T 3u <9T 1 - k 2 


/ uu\ 


LiAV =I-l—3- 




(A 8) 
(A 9) 

(A 10) 
(All) 

(A 12) 
(A 13) 
(A 14) 
(A 15) 
(A 16) 


(A 17) 
(A 18) 
(A 19) 
(A 20) 
(A21) 


(A 22) 


(p t duiiu 

In the case of isotropic background distributions fb(u), the potentials become functions of u only, 
and the friction vector and diffusion tensor can be simplified into 


~ , -u r 2 <9T 2 \ p 

K ab = ~ 47T— l ab y —- T^—\ ~ = ~ v l,abV, 


m b 


du c 2 du j p 


(A 23) 


n , ,,, 2y' ,:,T. Sy= 2T, 8 

D,„, = - i*T»y 14-,, - —^| - 


■ ^T ab y 


PP 


I^ + It 3 

u du c 2 


PP 


±^i + ±- 

uc 2 du c 4 


PP 

P 2 

I- 


PP 

rP 


= D, ab ^ + D, ab 1-^ 


(A 24) 
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The guiding-center transformation of the Fokker-Planck operator presented in (Brizard 2004) 
gave explicit expressions for the guiding-center friction and diffusion coefficients, \K“ b gc / and 

yDabgc}’ ' n t ^ le case °f isotropic background distributions and a non-relativistic collision ker¬ 
nel. Generalization of that work to relativistic collision kernel is straight-forward in the case of 
isotropic field-particle distributions because the forms of the particle phase-space friction and 
diffusion coefficients do not change. Only the expressions for v/ u /,, and i) u ,h are differ¬ 

ent but that will not affect the guiding-center transformation, as they are functions only of the 
guiding-center kinetic momentum. 
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